#################
# R????װ?????? #
#################
#install.packages(c("plyr", "dplyr", "xlsx", "multcomp", "ggplot2", "ggpubr", "RColorBrewer", "vegan"))
library(xlsx)                #??ȡд??excel
library(multcomp)            #???رȽ?
library(ggplot2)             #??ͼ
library(plyr)                #???ݿ?????
library(dplyr)               #???ݿ?????
library(ggpubr)              #??ͼ
library(RColorBrewer)        #??ͼ??ɫ
library(vegan)               #????????
library(ggrepel)             #???ǲ??ص???ǩ
library(ggThemeAssist)

############
# ???????? #
############
setwd("H://OneDrive//?о?//3-Сʵ??TU//5-????//TU?Ա??غ?????ֲ????Ӱ??//???ݷ???")  #???ݶ?ȡ???????ļ???
Sys.setlocale("LC_TIME", "English")
rm(list=ls())  #????ԭ?б?
F_name <- "TU_4N4I.xlsx"  #??ȡ??excel????????
A_start <- 5 #???????ڷ????ı???ʼ??
A_group <- c("Treat","PT") #????????ͳ?Ʒ????ı??? !!!-Ŀǰ????ֻ?ܷ???????-!!!
F_type <- 2 #??????ͼ???ͣ?1Ϊ????ͼ??2Ϊ????ͼ
F_x <- "Treat" #??????ͼ??x??
F_group <- "PT" #??????ͼ?ķ??飨ͼ????
F_ncol <- 0 #?????ϲ?ͼ??ÿ?е?Сͼ????0Ϊ????ͼÿ??ͼ??????????Ŀǰ??֧?ֺ?ͼ????
ylab <- list("Ht", "RL", "AB", "UB", "TB", "R/S", "LA", "LM", "SLA", "Pn", "WUE", "Chl", "Fv/Fm")
ytype <- list("Growth performance", "Growth performance", "Growth performance", "Growth performance", "Growth performance", "Growth performance", "Leaf traits", "Leaf traits", "Leaf traits", "Photosynthesis traits", "Photosynthesis traits", "Photosynthesis traits", "Photosynthesis traits")

windowsFonts(WF = windowsFont("Arial")) #???û?ͼ????

############
# ???ݴ??? #
############
data_t <- read.xlsx(F_name, sheetName = "T", encoding = "UTF-8")
data_o <- read.xlsx(F_name, sheetName = "data_o", encoding = "UTF-8")
group <- data_o[,A_group]
for (i in c(1:length(A_group))) {
  group[,i] <- factor(group[,i], levels = unique(group[,i]))     #?????????ӻ?
}
index_a <- colnames(data_o)[A_start:length(colnames(data_o))] #????ָ??#????ָ??

data_a <- list()
for (i in c(1:length(index_a))) {
  data_a[[i]] <- cbind(group, data_o[,index_a[i]])
  names(data_a[[i]]) <- c(A_group, "v")
  data_a[[i]]$v <- as.numeric(data_a[[i]]$v)
}


#????ͳ???㷨
result_F <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                     conf.interval=.95, .drop=TRUE) {
  library(plyr)
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=T) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2 (xx[[col]], na.rm=T),
                     mean = mean    (xx[[col]], na.rm=T),
                     sd   = sd      (xx[[col]], na.rm=T),
                     min  = min     (xx[[col]], na.rm=T),
                     max  = max     (xx[[col]], na.rm=T)
                   )
                 },
                 measurevar
  )
  datac$se <- datac$sd / sqrt(datac$N) 
  datac$down <- datac$mean - datac$se
  datac$up <- datac$mean + datac$se
  return(datac)
}



############
# ͳ?Ʒ??? #
############
###ͳ??????
result <- list()
for (i in c(1:length(index_a))) {
  result[[i]] <- result_F(data_a[[i]], measurevar = "v", groupvars = A_group)
  for (j in c(1:length(A_group))) {
    result[[i]][,j] <- as.numeric(factor(result[[i]][,j], levels = unique(result[[i]][,j])))
  }
}


#?????Լ???
for (i in c(1:length(index_a))) {                            #??????ָ??ѭ???ַ??????????б??ţ?
  star <- c("ns", "*", "**", "***", "****")
  t_P <- c(1,  0.05, 0.01, 0.001, 0.0001)
  for (j in c(1:length(A_group))) {                          #?????????????أ?ѭ???ַ???
    A <- as.numeric(group[,j])                               #????????A
    B <- as.numeric(group[,-j])                              #????????B
    if (sum(!duplicated(group[,j]))==2){                     #?ж????????????أ???ˮƽ?Ƿ?Ϊ2????Ϊ2??Ϊt????
      result_t <- data.frame()
      for (k in c(1:sum(!duplicated(B)))) {                #??????????Bˮƽѭ??????
        x1 <- data_a[[i]][A==1 & B==k,]$v 
        x2 <- data_a[[i]][A==2 & B==k,]$v                  #????????
        t <- t.test(x1, x2)[[1]]
        p <- t.test(x1, x2)[[3]]                           #t????
        for (l in c(1:5)) {
          if(p < t_P[l]){
            s <- star[l]
          }
        }
        result_t <- rbind(result_t,data.frame(k,t,p,s))    #??????¼
      }
      names(result_t) <- c(A_group[-j], 
                           paste(A_group[j], "t", sep = "_"), 
                           paste(A_group[j], "p", sep = "_"), 
                           paste(A_group[j], "d", sep = "_"))
      result[[i]] <- merge(result[[i]], result_t, by = A_group[-j], all.x = TRUE)
      rm(x1, x2, t, p, result_t)                             #?ϲ???????ͳ??????list????ɾ???м???
    } else {                                               #???ж????????????أ???ˮƽ??Ϊ2??Ϊ????????
      result_a <- data.frame()
      for (k in c(1:sum(!duplicated(B)))) {                #??????????Bˮƽѭ??????
        data_a[[i]] <- dplyr::rename(data_a[[i]], "factor" = A_group[j])
        #??????????A????Ϊ"factor"
        ANOVA <- aov(v ~ factor, data = data_a[[i]][B==k,])
        #????????
        data_a[[i]] <- plyr::rename(data_a[[i]], c("factor" = A_group[j]))
        #?ָ?????????A????
        tukey <- glht(ANOVA, linfct=mcp(factor = "Tukey")) #???رȽ?????ĸ
        result_a <- rbind(result_a, 
                          data.frame(c(1:sum(!duplicated(A))),
                                     rep(k,sum(!duplicated(B))),
                                     rep(summary(ANOVA)[[1]][1,4], sum(!duplicated(group[,j]))),
                                     rep(summary(ANOVA)[[1]][1,5], sum(!duplicated(group[,j]))),
                                     cld(tukey, level = .1)[[10]]$Letters))
      }  
      names(result_a) <- c(A_group[j], A_group[-j], paste(A_group[j], "F", sep = "_"), paste(A_group[j], "p", sep = "_"), paste(A_group[j], "d", sep = "_"))                        #????????????
      result[[i]] <- merge(result[[i]], result_a, by = c(A_group[-j], A_group[j]), all.x = TRUE)    
      rm(ANOVA, tukey, result_a)                           #?ϲ???????ͳ??????list????ɾ???м???
    }
    rm(A, B)
  }
}


#?????????ӻ?
for (i in c(1:length(index_a))) {
  result[[i]][,1] <- factor(result[[i]][,1], levels = unique(result[[i]][,1]))
  result[[i]][,2] <- factor(result[[i]][,2], levels = unique(result[[i]][,2]))
}



###????????????ֵ
data_re <- data.frame()
for (i in c(1:13)) {
  re <- cbind(data_a[[i]], ylab[i])
  names(re) <- c("Treat", "PT" ,"v", "Ind")
  data_re <- rbind(data_re, re)
}
for (i in c(1:length(data_re[,1]))) {
  data_re$v_re[i] <- (data_re$v[i] - result[[which(ylab == data_re$Ind[i])]]$mean[which(c("N", "I") == data_re$PT[i])])/result[[which(ylab == data_re$Ind[i])]]$mean[which(c("N", "I") == data_re$PT[i])]
}

data_re <- data_re[data_re$Treat != "T-U-",]
data_re$Ind <- factor(data_re$Ind, levels = rev(unique(data_re$Ind)))
result_re <- result_F(data_re, measurevar = "v_re", groupvars = c("Treat", "PT", "Ind"))

for (i in c(1:length(result_re[,1]))) {
  star <- c("ns", "*", "**", "***", "****")
  t_P <- c(1,  0.05, 0.01, 0.001, 0.0001)
  result_re[i,12] <- t.test(data_re$v_re[data_re$Treat == result_re$Treat[i] & data_re$Ind == result_re$Ind[i] & data_re$PT == "N"], data_re$v_re[data_re$Treat == result_re$Treat[i] & data_re$Ind == result_re$Ind[i] & data_re$PT == "I"])[[1]]
  result_re[i,13] <- t.test(data_re$v_re[data_re$Treat == result_re$Treat[i] & data_re$Ind == result_re$Ind[i] & data_re$PT == "N"], data_re$v_re[data_re$Treat == result_re$Treat[i] & data_re$Ind == result_re$Ind[i] & data_re$PT == "I"])[[3]]
  for (l in c(1:5)) {
    if(result_re[i,13] < t_P[l]){
      result_re[i,14] <- star[l]
    }
  }
  rm(star, t_P)
}
result_re <- plyr::rename(result_re, c("V12" = "t", "V13" = "P", "V14" = "s"))
result_re$Treat <- as.vector(result_re$Treat)
result_re$Treat[which(result_re$Treat == "T+U-")] <- "Temperature"
result_re$Treat[which(result_re$Treat == "T-U+")] <- "UV"
result_re$Treat[which(result_re$Treat == "T+U+")] <- "Temperature ?? UV"
result_re$Treat <- factor(result_re$Treat, levels = unique(result_re$Treat))
for (i in c(1:length(ylab))) {
  result_re$Ind_type[result_re$Ind == ylab[[i]]] <- ytype[[i]]
}
result_re$Ind_type <- factor(result_re$Ind_type, levels = rev(unique(result_re$Ind_type)))

########
# ??ͼ #
########
#????Ч??
p_T <- ggplot(data_t, aes(x=Date, y=Temperature, group=Treat, color=Treat)) +
  geom_line(size=1.2)+
  scale_x_date(date_breaks="10 day",date_labels="%m/%d")+
  labs(y= "Temperature (??)")+
  scale_y_continuous(limits = c(24, 39),
                     breaks = c(24, 29, 34, 39))+
  theme_bw()+
  theme(legend.position = "top",
        panel.grid=element_blank(),
        axis.title = element_text(size = 12, family = "WF"),
        axis.text =  element_text(size = 10, family = "WF"),
        strip.text = element_text(size = 10, family = "WF"),
        legend.text = element_text(size = 10, family = "WF"))+
  scale_color_manual(values=c("#4D93CA", "#E4636B"))

p <- list()
for (i in c(1:length(index_a))) {
  data_a[[i]] <- dplyr::rename(data_a[[i]], "x" = F_x, "g" = F_group)
  data_a[[i]]$x <- as.numeric(data_a[[i]]$x)
  result[[i]] <- dplyr::rename(result[[i]], "x" = F_x, "g" = F_group,
                               "text_x" = paste(F_x, "d", sep = "_"), 
                               "text_g" = paste(F_group, "d", sep = "_")) #????ͼ??x??????????׼??????
  
  p[[i]] <- ggboxplot(data_a[[i]], x="x", y="v", fill = "g")+
    stat_compare_means(method = "t.test", aes(group=g), label = "p.signif", size = 6)+
    geom_text(data = result[[i]], aes(label = text_x, color = g, x = x, y = -Inf), size = 6, 
              position = position_dodge(width=0.75), show.legend = FALSE, vjust = -0.2)+
    scale_fill_manual(values = c("#5B9BD5","#DE6757"),
                      breaks=c("N", "I"),
                      labels=c("A. argyi", "S. canadensis"))+
    scale_color_manual(values = c("#5B9BD5","#DE6757"))+
    scale_x_discrete(breaks = unique(result[[i]]$x),
                     labels = unique(group[,F_x]))+
    labs(x= "",y = ylab[i])+
    scale_y_continuous(expand = expansion(mult = c(0.2, 0.2)))+
    theme_bw()+
    theme(legend.position = "top",
          panel.grid = element_blank(),
          legend.title = element_blank(),
          axis.title = element_text(size = 16, family = "WF"),
          axis.text =  element_text(size = 14, family = "WF"),
          strip.text = element_text(size = 12, family = "WF"),
          legend.text = element_text(size = 14, family = "WF", face = "italic"))
  
  data_a[[i]] <- plyr::rename(data_a[[i]], c("x" = F_x, "g" = F_group))
  result[[i]] <- plyr::rename(result[[i]], c("x" = F_x, "g" = F_group,
                                             "text_x" = paste(F_x, "d", sep = "_"),
                                             "text_g" = paste(F_group, "d", sep = "_")))                    #?ָ?#????ͼ??x????????????
}




###
p_re <- ggplot(result_re)+
  geom_pointrange(aes(x = Ind, y = mean, ymax = up, ymin = down, color = PT, shape = PT), size=1.2, position = position_dodge(width=0.75))+
  geom_text(aes(label = s, x = Ind, y = -Inf), size = 8, show.legend = FALSE,hjust = -0.3)+
  labs(x= "",y = "")+
  scale_y_continuous(expand = expansion(mult = c(0.3,0.05)))+
  scale_color_manual(values = c("#5B9BD5","#DE6757"),
                     breaks=c("N", "I"),
                     labels=c("A. argyi", "S. canadensis"))+
  geom_hline(aes(yintercept = 0),linetype = 2)+
  theme_bw()+
  facet_grid(Ind_type ~ Treat, scales = "free_y", space = "free_y")+
  theme(legend.position="top",
        legend.title = element_blank(),
        panel.grid=element_blank(),
        axis.title = element_text(size = 20, family = "WF"),
        axis.text =  element_text(size = 18, family = "WF"),
        strip.text = element_text(size = 18, family = "WF"),
        legend.text = element_text(size = 18, family = "WF", face = "italic"))+
  guides(shape = "none")+
  coord_flip()



###????????
V_ex <- data_o[,c(1:4)]                              #??ȡ???ͱ?
V_re <- data_o[,c(5:length(data_o))]                 #??ȡ??Ӧ??
pt <- c("N", "I")
for (i in c(1:2)) {
  x <- V_ex[V_ex$PT == pt[i],][,c(2,3)]
  y <- V_re[V_ex$PT == pt[i],]
  result_rda <- summary(rda(y ~., x, scale=TRUE))
  sp <- as.data.frame(result_rda$species[,1:2])*5
  row.names(sp) <- ylab
  ###??ȡ??Ӧ?????꣬????5??ʹͼ?ۣ???Ӱ??????
  st <- cbind(V_ex[V_ex$PT == pt[i], 1], as.data.frame(result_rda$sites[,1:2]))
  names(st) <- c("Treat","RDA1","RDA2")
  ###??ȡ???????꣬????ģʽ???ɸ????Լ?????̽??????ѡһ????
  yz <- as.data.frame(result_rda$biplot[,1:2])*4
  ###??ȡ???ͱ?????
  xlab_RDA <- paste0("RDA1 (", round(result_rda[["cont"]][["importance"]][2,1]*100, 2), "%)")
  ylab_RDA <- paste0("RDA2 (", round(result_rda[["cont"]][["importance"]][2,2]*100, 2), "%)")
    p[[length(p)+1]] <- ggplot() +
      geom_point(data = st, aes(x = RDA1, y = RDA2, colour = Treat, shape = Treat), 
                 size = 3, stroke = 2)+
      stat_ellipse(data = st, aes(x = RDA1, y = RDA2, fill = Treat),
                   type = "t", geom = "polygon",alpha = 0.1,color = NA)+
      scale_shape_manual(values = c(0, 1, 15, 16))+
      scale_colour_manual(values = brewer.pal(4,"Set2"))+
      scale_fill_manual(values = brewer.pal(4,"Set2"))+
      geom_segment(data = sp,aes(x = 0, y = 0, xend = RDA1, yend = RDA2),
                   arrow = arrow(angle=22.5,length = unit(0.35,"cm"), type = "closed"),
                   linetype = 1, size = 0.8, colour = brewer.pal(8,"Set1")[1])+ 
      geom_text_repel(data = sp, aes(RDA1, RDA2, label = row.names(sp)),
                      colour = brewer.pal(8,"Set1")[1], position = "identity")+
      geom_segment(data = yz,aes(x = 0, y = 0, xend = RDA1, yend = RDA2),
                   arrow = arrow(angle=22.5,length = unit(0.35,"cm"), type = "closed"),
                   linetype = 1, size = 1.5, colour = brewer.pal(8,"Set1")[2])+  
      geom_text_repel(data = yz, aes(RDA1, RDA2, label = row.names(yz)),
                      colour = brewer.pal(8,"Set1")[2], position = "identity")+
      geom_text(aes(x = c(5,5,-5,-5), y = c(5,-5,-5,5), label = c("(I)", "(II)", "(III)", "(IV)")), show.legend = F, size = 6, family = "WF")+
      labs(x = xlab_RDA, y= ylab_RDA)+
      geom_hline(yintercept = 0,linetype = 3,size = 1)+ 
      geom_vline(xintercept = 0,linetype = 3,size = 1)+
      scale_x_continuous(limits = c(-5, 5))+
      scale_y_continuous(limits = c(-5, 5))+
      theme_bw()+
      theme(panel.grid = element_blank(),
            legend.title = element_blank(),
            axis.title = element_text(size = 16, family = "WF"),
            axis.text =  element_text(size = 14, family = "WF"),
            strip.text = element_text(size = 14, family = "WF"),
            legend.text = element_text(size = 14, family = "WF"))
  result[[length(result)+1]] <- result_rda[["species"]]
  result[[length(result)+1]] <- result_rda[["cont"]][["importance"]]
}



############
# ???????? #
############
sheetName <- c(index_a, "Native_sp", "Native_Im", "Invaisve_sp", "Invaisve_Im")
##??????
write.xlsx(result[[1]],"????????.xlsx", sheetName = sheetName[1])
if (length(result) > 1) {
  for (i in c(2:length(result))) {
    write.xlsx(result[[i]],"????????.xlsx", sheetName = sheetName[i], append = TRUE)
  }
}
write.xlsx(result_re,"????????.xlsx", sheetName = "re", append = TRUE)

##????ͼ
tiff(file = "000.tiff", width = 2000, height = 800, units = "px", res = 300)
print(p_T)
dev.off()

tiff(file = "001.tiff", width = 2200, height = 2500, units = "px", res = 300)
print(ggpubr::ggarrange(p[[1]],p[[3]],p[[2]],p[[4]],p[[6]],p[[5]], nrow = 3, ncol = 2, 
                        labels = LETTERS[1:6], vjust = 0.5,
                        font.label = list(color = 'red'), common.legend = TRUE, align = "hv"))
dev.off()

tiff(file = "002.tiff", width = 1100, height = 2500, units = "px", res = 300)
print(ggpubr::ggarrange(p[[7]],p[[8]],p[[9]], nrow = 3, ncol = 1, 
                        labels = LETTERS[1:6], vjust = 0.5,
                        font.label = list(color = 'red'), common.legend = TRUE, align = "hv"))
dev.off()

tiff(file = "003.tiff", width = 2200, height = 1600, units = "px", res = 300)
print(ggpubr::ggarrange(p[[10]],p[[11]],p[[12]],p[[13]], nrow = 2, ncol = 2, 
                        labels = LETTERS[1:6], vjust = 0.5,
                        font.label = list(color = 'red'), common.legend = TRUE, align = "hv"))
dev.off()

tiff(file = "004.tiff", width = 3000, height = 1600, units = "px", res = 300)
print(ggpubr::ggarrange(p[[14]],p[[15]], nrow = 1, ncol = 2, 
                        labels = c("A. argyi", "S. canadensis"), vjust = 0.2,
                        font.label = list(color = 'red', face = "italic"), 
                        common.legend = TRUE, align = "hv"))
dev.off()

tiff(file = "S1.tiff", width = 3100, height = 3000, units = "px", res = 300)
print(p_re)
dev.off()

